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We present a numerical scheme for solving the time-independent nonlinear Gross-Pitaevskii equation 
in two dimensions describing the Bose-Einstein condensate of trapped interacting neutral atoms 
at zero temperature. The trap potential is taken to be of the harmonic-oscillator type and the 
interaction both attractive and repulsive. The Gross-Pitaevskii equation is numerically integrated 
consistent with the correct boundary conditions at the origin and in the asymptotic region. Rapid 
convergence is obtained in all cases studied. In the attractive case there is a limit to the maximum 
number of atoms in the condensate. 
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Recently, there have been experiments [y of Bose- 
Einstein condensation in dilute bosonic atoms (alkali and 
hydrogen atoms) employing magnetic traps at ultra-low 
temperatures. These experiments have intensified theo- 
retical investigations on various aspects of the condensate 
. The condensate can consist of few thousand to mil- 
lions of atoms confined by the trap potential. The prop- 
erties of the condensate at zero temperature are usually 
described by the nonlinear time-independent mean-field 
Gross-Pitaevskii (GP) equation which properly in- 
corporates the trap potential as well as the interaction 
among the atoms. The effect of the interaction leads to 
a nonlinear term in the GP equation which complicates 
its solution procedure. 

There have been a series of recent studies which deal 
with the numerical solution of the three-dimensional GP 
equation |^-|^. These works have emphasized the seri- 
ous difficulties in obtaining numerical convergence of the 
solution. There has been no such systematic study on 
the numerical solution of the GP equation in two dimen- 
sions. Although, there has been no experiments on Bose- 
Einstein condensation of two-dimensional systems, this 
is a problem of great interest. A system of ideal Bose 
gas in two dimensions does not undergo condensation 
at a finite temperature |9|. However, condensation can 
take place under the action of a trap potential In 
order to achieve two-dimensional Bose-Einstein conden- 
sation in real three-dimensional traps, one should choose 
the frequency ojxy in the x — y plane negligibly small 
compared to that in the z direction ujz [0. Also, there 
has been consideration of Bose-Einstein condensation in 
low-dimensional systems for particles confined by gravi- 
tational field or by a rotational container Possible 
experimental configurations for Bose-Einstein condensa- 
tion in spin-polarized hydrogen in two dimensions are 
currently being discussed Because of these inter- 

ests, here we perform a critical study of the numerical 
solution of the time-independent GP equation in two di- 
mensions for an interacting Bose gas under the action of a 
harmonic-oscillator-type trap potential. The interatomic 
interaction is taken to be both attractive and repulsive 



in nature. 

In steady state at zero temperature the condensate 
wave function is described by the following effective non- 
linear Schrodinger-type equation known as the Gross- 
Pitaevskii equation for condensed neutral bosons in 
a harmonic trap: 
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-t- -mcj^r^ -I- o^'^(r) - a 



^'(r)=0. (0.1) 



Here 4'(r) is the condensate wave function at position r, 
m the mass of a single bosonic atom, muPr^ /2 the attrac- 
tive harmonic-oscillator trap potential, ld the oscillator 
frequency, /i the chemical potential and g the strength 
of interatomic interaction. A positive g correspond to 
a repulsive interaction and a negative g to an attractive 
interaction. 

Here we shall be interested i n th e spherically symmet- 
ric solution ^'(r) = ?/'(r) to eq. (0.1) which can be written 
as 



Tt? I d d 1 2 2 ,2/ N 



i/'(r) = 0. (0.2) 



The ground state of the condensate appears in such a 
spherically symmetric st ate. As is Ref. it is con- 
venient to express eq. (0.2) in terms of dimensionless 
variables defined hy x ^ r/a, where a = \/h/ (muj), 
a — ii/{huj), ip^x) — ay'2 mg'ip( r) /h. In terms of these 
dimensionless variables eq. (|0.2|) can be written as 



—X— — ha; + c^p (x) — 2a 

X dx dx 



tpix) = 0. (0.3) 



where c = ±1 carries the sign of g, c = 1 corresponds 
to a repulsive interaction and c = — 1 corresponds to an 
attractive interaction. 

The normalization of the wave function is given by 



N = 2t: 



drr^p'^{r), 



(0.4) 
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where N is the total number of atoms in the condensate. 
In terms of the dimensionless variables defined above, 
this normalization condition becomes 



xdxtp'^{x) — n = j]N, 



(0.5) 



where as in Ref. we have introduced a dimensionless 
coupling T] = mg/(TTh^) of interatomic interactions and 
a reduced number of particles n. From a study of the 
temperature dependence of the chemical potential for a 
system of interacting bosons in two dimensions it was 
concluded in Ref. that a Bose-Einstein-condensate-like 
behavior is obtained in this system for values of 77 typi- 
cally smaller than 0.001. Hence for a qualitative estimate 
of the number of condensed bosons in this calculation, we 

shall consider = 0.0001. 

Instead of solvin g eq s. (|0.2|) and (0.4), w e sh all be 
working with eqs. ( p.3[) and ( p.5[ ). Equation (X3) is in- 



dependent of all parameters of the problem, such as, to, 
uj, g, and N. The rele vant parameters appear in the nor- 
malization condition (0.5). The constant n in eq. ( |0.5| ) 
is the reduced number of atoms for the system and is 
related to the real number of atoms N. 

Another interesting property of the condensate wave 
function is its mean-square radius defined by 



(0.6) 



In terms of the dimensionless variables defined above we 
have 



(.T^) = — I x'^^p^{x)xdx. 



(0.7) 



To solve eq. ( p.BD numerically, first we study the 
asymptotic behavior of its solutions. Since for a suf- 
ficiently large x, ip{x) must vanish asymptotically, the 
nonlinear term p rop ortional to ^l'^{x) can eventually be 
neglected in eq. (0^) in the asymptotic region. Thus for 
large x, this equation can be approximated by 



1 d d 2 
X dx dx 



2a 



Tp{x) = 0. 



(0.8) 



If this equation would be valid for all x this would 
be the equation for the two-dimensional oscillator in 
the spherically symmetric state permitting solution for 
a = 1,3,5,... etc. However, in the present problem eq. 
( |0.8|) is valid only in the asymptotic region. Considering 
eq. ( p.SD as a mathematical equation valid for all a, the 
asymptotic form of the physically acceptable solution is 
given by 



lim ip{x) — Nc exp 



— + (a — 1) In a; 



(0.9) 



where Nc is a normalization constant. The derivative 
of the wave function in the asymptotic region can be 



obtained from eq. (0.9) and one obtains the following 
asymptotic form for the log-derivative of the wave func- 
tion 



lim 



tp{x) 



a-1 



(0.10) 



which i s in dependent of the normalization constant Nc 
of eq. (3.9) and where the prime denotes derivative with 

respect to x. 

Next we consider eq. (0.3) as x — > 0. The nonlinear 
term approaches a constant in this limit because of the 
regularity of the wave function at x = 0. Then one has 
the following usual conditions 



■0(0) = constant, 



^'(0) = 0, 



(0.11) 



as in the case of the two-dimensional harmonic oscillator 
problem. 



Equation (0.3) is integrated numerically for a given a 
by the four-point Runge-Kutta rule starting at the origin 
{x = 0) with the initial boundary condition (0.11) with 
a trial V'(O). The numerical integration is performed in 
steps of dx = A, where A is typically taken to be 0.0001. 
The integration is pro paga ted to x = Xmax, where the 
asymptotic condition ( |0.10| ) is valid. The agreement be- 
tween the numerically calculated log-der ivativ e of the 
wave function and the theoretical result (3.10) was en- 
forced to five significant figures. The maximum value of 
X, up to which we needed to integrate (3.3) numerically 
for obtaining this precision, is Xmax — 5. If for a trial 
^(0), the agreement of the log-derivative can not be ob- 
tained, a new value of '0(0) is to be chosen. The proper 
choice of ip{0) was implemented by the secant method. 
Even with this method, sometimes it is difficult to ob- 
tain the proper value of 0(0) for a given a. Unless the 
initial guess is "right" and one is sufhciently near the de- 
sired solution the method could fail and lead to either 
the trivial solution 0'(x) = or an exponentially diver- 
gent nonnormalizable solution in the asymptotic region. 
However, there is one guideline which is of help in finding 
the proper value of 0(0). The convergence of the numer- 
ical log-derivative to the theoretical result (3.10) with a 
variation of V'(O) for a fixed a is monotonic in nature 
provided that one is near the exact value. If for two trial 
values of 0(0) near the exact value one obtains a positive 
and a negative error for the log-derivative, respectively, 
the exact 0(0) lies in between these two trial values. 

For both attractive (c = —1) and repulsive (c = 1) 
interatomic interactions, in ad diti on to the ground state 
solution with no nodes, eq. (13) also permits radially 
excited states with nodes in the wave function. In the 
absence of the nonlinear term, the spherically-symmetric 
discrete ground-state solu tion of the two-dimensional os- 
cillator given by eq. (3.3) is obtained for a = 1. Radially 
excited oscillator states appear for a = 3, 5, ... etc. In the 
presence of the nonlinearity, for attractive interatomic 
interaction (c = —1), the solutions of the GP equation 
for the ground state appear for values chemical potential 
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a < 1. For a repulsive interatomic interaction (c = 1), 
t hes e solutions appear for a > 1. The solutions of eq. 
( |0.3[ ) for the radially excited state with one node for the 
attractive (c — —1) and repulsive (c = 1) cases appear for 
values of chemical potential a smaller and larger than the 
energy of the first excited state {a — 3) of the harmonic 
oscillator problem, respectiwrfy. 



ture of the wave function for these two cases are quite 
different. 

Table 1: Pa ram eters for the numerical solution of the 
GP equation (13) for c = ±1 for the ground state wave 
function. The first four columns refer to the attractive 
interaction c = — 1 and the last four columns refer to the 
repulsive interaction c = 1. 
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Fig. 1. Ground-state condensate wave function tp{x) 
versus x for (a) attractive and (b) repulsive interparticle 
interactions. The parameters for these cases are given in 
Table I. The curves appear in the same order as in Tabic 
I. The lowermost curve corresponds to the first row in 
Table I. 

First we consider the ground-state solution of eq. ( |0.3| ) 
for different a in the cases of both attractive and repulsive 
interactions. The relevant parameters for these solutions 
(values of the wave-function at the origin ip{0), reduced 
number n, and mean-square radii (x^)) are listed in ta- 
ble 1. The wave functions for different values of a for 
the attractive and repulsive interparticle interactions for 
the cases shown in table 1 are exhibited in figures 1(a) 
and 1(b), respectively, where we plot tpix) versus x. The 
curves in figures 1(a) and 1(b) appear in the same order 
as the rows in table 1 and it is easy to identify the cor- 
responding values of a from the values of ipW oi each 
curve. From figures 1(a) and (b) we find that the na- 
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Fig. 2. The reduced number of particles n versus 
density p = n/{x'^) of the ground-state condensate for 
the attractive interaction. 

Table 2: Pa ram eters for the numerical solution of the 
GP equation ( |0.3|) for c — ±1 for the wave function of 
the first excited state. The first four columns refer to the 
attractive interaction c — —1 and the last four columns 
refer to the repulsive interaction c = 1. 
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For the attractive interparticle interaction, the wave 
function is more sharply peaked at a; = than in the case 
of the repulsive interparticle interaction. Consequently, 
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in the attractive case one has smaUer values of the re- 
duced number n and mean square radii (x^) of the wave 
function. For the attractive case we find from table 1 
that with a reduction of the chemical potential a the re- 
duced number n increases slowly and the mean square 



radius 



decreases rapidly, so that the density of the 



condensate p = n/{x'^) tends to diverge as n tends to a 
maximum value nmax- This means that there is a max- 
imum number of particles in the condensate in this case. 
This peculiar behavior in the attractive case is demon- 
strated in figure 2 where we plot the reduced number n 
of the condensate in the ground state versus the density 
p. From figure 2 we find that with the reduction of the 
chemical potential, the reduced number of particles n in 
the condensate attains a saturation value. Numerically, 
we find this value to be 



"max = ?7^max ~ 1-85. 



(0.12) 



However, during this process the mean square radius con- 
tinues to decrease thus leading to a divergence of the den- 
sity of the condensate. For n > riniax, there is no stable 
solution of the GP equation. There is no such limit on 
n in the repulsive case. In that case with the increase of 
the chemical potential a the condensate increases in size 
as the number of particles in the condensate increases. 
These behaviors of the Bose-Einstein condensate in two 
dimensions were also noted in three dimensions ^1 13 1. 

Next we consider solutions to eq. (0.2) with radial 
excitation. The present numerical method is equally ap- 
plicable to the ground and all radially excited states of 
the condensate. In this work we shall consider excited 
solutions with only one node for both attractive and re- 
pulsive interparticle interactions, which we discuss next. 
In table 2 we exhibit the parameters of some of such solu- 
tions. In figures 3(a) and 3(b) we plot the respective wave 
functions ipi^) versus x for the attractive and repulsive 
interactions. We again find that in the case of attractive 
interaction the wave function is narrowly peaked near 
X = 0, and for repulsive interaction it is more extended 
in space to larger values of x. Consequently, the reduced 



number n and mean square radii 



are larger in the 



case of repulsive interaction. We again find from table 2 
that in the attractive case, with a reduction of the chem- 
ical potential a, the mean square radius decreases as the 
reduced number n increases. It is expected that there 
should be a saturation on n in this case also. However, 
we did not try to find this saturation numerically, which 
should occur for values of a much smaller than those pre- 
sented in table 2 leading to a much larger value of nniax 
than in eq. ( |d1^ ). 

Although, we have considered the problem in a sys- 
tem of dimensionless units, it is interesting to see what 
our results correspond to in actual units. If we consider 
the typical value 0.0001 of rj favorable to the formation 
of B ose-Einstein condensate Q as commented after eq. 
( p.SD , we can calculate the actual numer of atoms N. 
The number of atoms N in the condensate for all cases 



reported in tables 1 and 2 is maximum for the ground 
state of the repulsive interatomic interaction for a — -i. 
This number in this case is = 147, 600. If we consider 
a trap frequency such that y/h/Jmuj) — 10,000 A 
then the root-mean-square radius for the above conden- 
sate with 147,600 atoms is 16,700 A. In the attr active 
case the maximum numer of particles given byeq. dol^ ) 
becomes A^max = 18, 500. Both the size and numbers of 
particles seem to be very reasonable for the condensate 
i|. 

In this work we have investigated the numerical so- 
lution of the Gross-Pitaevskii equation (0.1) for Bose- 
Einstein condensation in two dimensions under the action 
of a harmonic oscillator trap potential for bosonic atoms 
interacting via both attractive and repulsive interparticle 
interactions. In both cases we considered the wave func- 
tion for the ground state and radially excited state with 
one node. We expressed the GP equation in dimension- 
less units independent of all parameters, such as, atomic 
mass, harmonic oscillator frequency, number of atoms in 
the condensate, and strength of atomic interaction. The 
relev ant parameters appear in the normalization condi- 
tion (0.5) of the wave hmction. We derive the boundary 
conditions (0.10) and (|0.11|) of the solution of the di- 



mensionless GP equation ( |0.3| ), which is integrated from 
the origin outwards in steps of 0.0001 by the four-point 
Runge-Kutta rule consistent with the boundary condi- 
tion. At a particular value of the chemical potential, the 
correct solution is obtained after a proper guess of the 
boundary condition of the wave function at the origin. 
From a initial trial value of the wave function at the ori- 
gin, Newton-Raphson method is used to obtain the cor- 
rect wave function after a matching with the boundary 
condition in the asymptotic region. In both the ground 
and excited states it is found that the wave function is 
sharply peaked near the origin for attractive interatomic 
interaction. For a repulsive interatomic interaction the 
wave function extends over a larger region of space. In 
the case of an attractive potential, the mean square ra- 
dius decreases with an increase of the number of particles 
in the condensate. Consequently, a stable solution of the 
GP equation can be obtained for a maximum number 
of particles in the condensate. For the ground state the 
maxim um reduced number of particles is given by eq. 
( |0.12 ). For the repulsive case there is no such limit on 



the number of particles in the condensate. 
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